* Set Directory
clear
set more off

cd "$path"
global working_data "$path/data"
global out_files "$path/output"

********************************************************************************
* Table 1: Number of observations for unemployed men by country and Early Retirement Age (also Appendix Table C.3 for women)
********************************************************************************	
clear
use $working_data/SHARE_JOLE_replication.dta, clear

keep if inrange(ragey, 50, 70) // only keep individuals aged 50-70

*ERA in years
gen group_eray = round(era_months/12, 1.0)

*Define artificial never treated group:
gen group_eray_never = group_eray
replace group_eray_never = 0 if group_eray==63 | group_eray==64

label define group_eray_label 55 "55" 56 "56" 57 "57" 58 "58" 59 "59" 60 "60" 61 "61" 62 "62" 63 "63" 64 "64" 
label value group_eray group_eray_label

*Men (Table 1)
jwdid depression if ragey<63 & ragender==1 & c_lms_new==3, ivar(nmergeid) tvar(ragey) gvar(group_eray_never) never fevar(riwm)

latab isocountry group_eray if e(sample), tf($out_files/crosstable_country_ERA_men) replace

*Women (Appendix Table C.3)
jwdid depression if ragey<63 & ragender==2 & c_lms_new==3, ivar(nmergeid) tvar(ragey) gvar(group_eray_never) never fevar(riwm)

latab isocountry group_eray if e(sample), tf($out_files/crosstable_country_ERA_women) replace

********************************************************************************
* Figure 1: Mental well-being between ages 50-70 (also Appendix Figures C.1 and C.3)
********************************************************************************
*Note: the following code creates the panels for Figure 1 as well as Appendix Figures C.1 and C.3. In addition the code generates the equivalent of Appendix Figure C.1 for women, which is not included in the paper's results

clear
use $working_data/SHARE_JOLE_replication.dta, clear

gen d_unemp = c_lms_new==3
gen d_disab = c_lms_new==4

*********************
* Depression score
*********************
*Loop over gender
levels ragender, local(ragender)
local forgender: value label ragender
foreach x of local ragender{
local gender_label: label `forgender' `x'

	matrix depression_all_`gender_label' = J(21, 5, .)

	*EURO-D scores for all countries combined
	forvalues i=50(1)70{
		
			matrix depression_all_`gender_label'[`i'-49,1]=`i'
			
			*Employed/Retired
			sum depression if ragender==`x' & ragey==`i' & (c_lms_new==1 | c_lms_new==2)
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix depression_all_`gender_label'[`i'-49,2]=mental_health
			
			*Unemployed
			sum depression if ragender==`x' & ragey==`i' & c_lms_new==3
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix depression_all_`gender_label'[`i'-49,3]=mental_health
			
			*Disabled
			sum depression if ragender==`x' & ragey==`i' & c_lms_new==4
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix depression_all_`gender_label'[`i'-49,4]=mental_health

			*Full population
			sum depression if ragender==`x' & ragey==`i' 
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix depression_all_`gender_label'[`i'-49,5]=mental_health

		}
			
	svmat depression_all_`gender_label'

	*Figure: Employed/Retired vs. Unemployed
	reg depression i.d_unemp d_unemp#c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	if (c_lms_new==1 | c_lms_new==2 | c_lms_new==3) /// 
	& ragender==`x' & ragey>=50 & ragey<=70 , vce(cluster nmergeid)

			mat b = e(b)
			mat c = e(V)
			
			scalar v=b[1,3]
			local coef1=v
			local coef1 : di %4.3f `coef1'
			
			scalar v_var=c[3,3]
			local st_err1=sqrt(v_var)
			local st_err1 : di %4.3f `st_err1'
			
			scalar z=b[1,4]
			local coef2=z
			local coef2 : di %4.3f `coef2'
			
			scalar z_var=c[4,4]
			local st_err2=sqrt(z_var)
			local st_err2 : di %4.3f `st_err2'

	margins, at(ragey = (50(1)70) d_unemp = (0 1))

	marginsplot,  recastci(rarea) ci1opts(fintensity(30)) recast(line) xlabel(50(2)70) ylabel(.1(.05).35, angle(0)) ///
	title(Slope E/R: `coef1' (`st_err1'); Slope U: `coef2' (`st_err2'), color(black) size(medium) ) /// 
	addplot(scatter depression_all_`gender_label'3 depression_all_`gender_label'1, mcolor(black) msize(2) msymbol(X) xlabel(50(2)70) || scatter depression_all_`gender_label'2 depression_all_`gender_label'1, mcolor(black) msize(2) msymbol(t) xlabel(50(2)70) ylabel(.1(.05).35, angle(0))) ///
	xtitle("Age") ytitle("Fraction depressed") /// 
	legend(off) ///
	name(conv_depression_`gender_label')
	
	graph display conv_depression_`gender_label', xsize(1.5) ysize(1) scale(1.3)
	
	graph export "$out_files/conv_depression_`gender_label'.pdf", replace
	
	*Figure: Employed/Retired vs. Disabled
	reg depression i.d_disab d_disab#c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	if (c_lms_new==1 | c_lms_new==2 | c_lms_new==4) /// 
	& ragender==`x' & ragey>=50 & ragey<=70 , vce(cluster nmergeid)

			mat b = e(b)
			mat c = e(V)
			
			scalar v=b[1,3]
			local coef1=v
			local coef1 : di %4.3f `coef1'
			
			scalar v_var=c[3,3]
			local st_err1=sqrt(v_var)
			local st_err1 : di %4.3f `st_err1'
			
			scalar z=b[1,4]
			local coef2=z
			local coef2 : di %4.3f `coef2'
			
			scalar z_var=c[4,4]
			local st_err2=sqrt(z_var)
			local st_err2 : di %4.3f `st_err2'

	margins, at(ragey = (50(1)70) d_disab = (0 1))

	marginsplot,  recastci(rarea) ci1opts(fintensity(30)) recast(line) xlabel(50(2)70) ylabel(1.5(.5)4, angle(0)) ///
	title(Slope E/R: `coef1' (`st_err1'); Slope D: `coef2' (`st_err2'), color(black) size(medium) ) /// 
	addplot(scatter depression_all_`gender_label'4 depression_all_`gender_label'1, mcolor(black) msize(2) msymbol(X) xlabel(50(2)70) || scatter depression_all_`gender_label'2 depression_all_`gender_label'1, mcolor(black) msize(2) msymbol(t) xlabel(50(2)70) ylabel(.1(.1).7, angle(0))) ///
	xtitle("Age") ytitle("Fraction depressed") /// 
	legend(off) ///
	name(convD_depression_`gender_label')
	
	graph display convD_depression_`gender_label', xsize(1.5) ysize(1) scale(1.3)
	
	graph export "$out_files/convD_depression_`gender_label'.pdf", replace
	
	*Figure: Full population
	reg depression c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry i.c_lms_new ///
	if ragender==`x' & ragey>=50 & ragey<=70 , vce(cluster nmergeid)

			mat b = e(b)
			mat c = e(V)
			
			scalar v=b[1,1]
			local coef1=v
			local coef1 : di %4.3f `coef1'
			
			scalar v_var=c[1,1]
			local st_err1=sqrt(v_var)
			local st_err1 : di %4.3f `st_err1'
			
	margins, at(ragey = (50(1)70))

	marginsplot,  recastci(rarea) ci1opts(fintensity(30)) recast(line) xlabel(50(2)70) ylabel(.1(.05).35, angle(0)) ///
	title(Slope: `coef1' (`st_err1'), color(black) size(medium) ) /// 
	addplot(scatter depression_all_`gender_label'5 depression_all_`gender_label'1, mcolor(black) msize(2) xlabel(50(2)70) ylabel(.1(.05).35, angle(0))) ///
	xtitle("Age") ytitle("Fraction depressed") /// 
	legend(off) ///
	name(pop_depression_`gender_label')
	
	graph display pop_depression_`gender_label', xsize(1.5) ysize(1) scale(1.3)
	
	graph export "$out_files/pop_depression_`gender_label'.pdf", replace	
}

*********************
* EURO-D score
*********************
*Loop over gender
levels ragender, local(ragender)
local forgender: value label ragender
foreach x of local ragender{
local gender_label: label `forgender' `x'

	matrix eurod_all_`gender_label' = J(21, 5, .)

	*EURO-D scores for all countries combined
	forvalues i=50(1)70{
		
			matrix eurod_all_`gender_label'[`i'-49,1]=`i'
			
			*Employed/Retired
			sum reurod if ragender==`x' & ragey==`i' & (c_lms_new==1 | c_lms_new==2)
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix eurod_all_`gender_label'[`i'-49,2]=mental_health
			
			*Unemployed
			sum reurod if ragender==`x' & ragey==`i' & c_lms_new==3
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix eurod_all_`gender_label'[`i'-49,3]=mental_health
			
			*Disabled
			sum reurod if ragender==`x' & ragey==`i' & c_lms_new==4
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix eurod_all_`gender_label'[`i'-49,4]=mental_health
			
			*Full population
			sum reurod if ragender==`x' & ragey==`i' 
			scalar total=r(N)
			if total<15 scalar mental_health=.
			if total>=15 scalar mental_health=r(mean)
			matrix eurod_all_`gender_label'[`i'-49,5]=mental_health

		}
			
	svmat eurod_all_`gender_label'

	*Figure: Employed/Retired vs. Unemployed
	reg reurod i.d_unemp d_unemp#c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	if (c_lms_new==1 | c_lms_new==2 | c_lms_new==3) /// 
	& ragender==`x' & ragey>=50 & ragey<=70 , vce(cluster nmergeid)

			mat b = e(b)
			mat c = e(V)
			
			scalar v=b[1,3]
			local coef1=v
			local coef1 : di %4.3f `coef1'
			
			scalar v_var=c[3,3]
			local st_err1=sqrt(v_var)
			local st_err1 : di %4.3f `st_err1'
			
			scalar z=b[1,4]
			local coef2=z
			local coef2 : di %4.3f `coef2'
			
			scalar z_var=c[4,4]
			local st_err2=sqrt(z_var)
			local st_err2 : di %4.3f `st_err2'

	margins, at(ragey = (50(1)70) d_unemp = (0 1))

	marginsplot,  recastci(rarea) ci1opts(fintensity(30)) recast(line) xlabel(50(2)70) ylabel(1.5(.5)4, angle(0)) ///
	title(Slope E/R: `coef1' (`st_err1'); Slope U: `coef2' (`st_err2'), color(black) size(medium) ) /// 
	addplot(scatter eurod_all_`gender_label'3 eurod_all_`gender_label'1, mcolor(black) msize(2) msymbol(X) xlabel(50(2)70) || scatter eurod_all_`gender_label'2 eurod_all_`gender_label'1, mcolor(black) msize(2) msymbol(t) xlabel(50(2)70) ylabel(1.5(.25)3, angle(0))) ///
	xtitle("Age") ytitle("EURO-D score") /// 
	legend(off) ///
	name(conv_eurod_`gender_label')
	
	graph display conv_eurod_`gender_label', xsize(1.5) ysize(1) scale(1.3)
	
	graph export "$out_files/conv_eurod_`gender_label'.pdf", replace
	
	*Figure: Employed/Retired vs. Disabled
	reg reurod i.d_disab d_disab#c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	if (c_lms_new==1 | c_lms_new==2 | c_lms_new==4) /// 
	& ragender==`x' & ragey>=50 & ragey<=70 , vce(cluster nmergeid)

			mat b = e(b)
			mat c = e(V)
			
			scalar v=b[1,3]
			local coef1=v
			local coef1 : di %4.3f `coef1'
			
			scalar v_var=c[3,3]
			local st_err1=sqrt(v_var)
			local st_err1 : di %4.3f `st_err1'
			
			scalar z=b[1,4]
			local coef2=z
			local coef2 : di %4.3f `coef2'
			
			scalar z_var=c[4,4]
			local st_err2=sqrt(z_var)
			local st_err2 : di %4.3f `st_err2'

	margins, at(ragey = (50(1)70) d_disab = (0 1))

	marginsplot,  recastci(rarea) ci1opts(fintensity(30)) recast(line) xlabel(50(2)70) ylabel(1.5(.5)4, angle(0)) ///
	title(Slope E/R: `coef1' (`st_err1'); Slope D: `coef2' (`st_err2'), color(black) size(medium) ) /// 
	addplot(scatter eurod_all_`gender_label'4 eurod_all_`gender_label'1, mcolor(black) msize(2) msymbol(X) xlabel(50(2)70) || scatter eurod_all_`gender_label'2 eurod_all_`gender_label'1, mcolor(black) msize(2) msymbol(t) xlabel(50(2)70) ylabel(1.5(1)5, angle(0))) ///
	xtitle("Age") ytitle("EURO-D score") /// 
	legend(off) ///
	name(convD_eurod_`gender_label')
	
	graph display convD_eurod_`gender_label', xsize(1.5) ysize(1) scale(1.3)
	
	graph export "$out_files/convD_eurod_`gender_label'.pdf", replace
	
	*Figure: Full population
	reg reurod c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry i.c_lms_new ///
	if ragender==`x' & ragey>=50 & ragey<=70 , vce(cluster nmergeid)

			mat b = e(b)
			mat c = e(V)
			
			scalar v=b[1,1]
			local coef1=v
			local coef1 : di %4.3f `coef1'
			
			scalar v_var=c[1,1]
			local st_err1=sqrt(v_var)
			local st_err1 : di %4.3f `st_err1'
	
	margins, at(ragey = (50(1)70))

	marginsplot,  recastci(rarea) ci1opts(fintensity(30)) recast(line) xlabel(50(2)70) ylabel(1.5(.5)4, angle(0)) ///
	title(Slope: `coef1' (`st_err1'), color(black) size(medium) ) /// 
	addplot(scatter eurod_all_`gender_label'5 eurod_all_`gender_label'1, mcolor(black) msize(2) xlabel(50(2)70) ylabel(1.5(.25)3, angle(0))) ///
	xtitle("Age") ytitle("EURO-D score") /// 
	legend(off) ///
	name(pop_eurod_`gender_label')
	
	graph display pop_eurod_`gender_label', xsize(1.5) ysize(1) scale(1.3)
	
	graph export "$out_files/pop_eurod_`gender_label'.pdf", replace	
	
}

********************************************************************************
* Appendix Table F.1: Estimated linear age coefficients for mental well-being of men between ages 50-70
********************************************************************************
clear all
use $working_data/SHARE_JOLE_replication.dta, clear

local file country_age_table.tex

gen d_unemp = c_lms_new==3
gen d_disab = c_lms_new==4

*Loop over gender
levels ragender, local(ragender)
local forgender: value label ragender
foreach x of local ragender{
local gender_label: label `forgender' `x'

*Depression
writeln $out_files/`file' "\multicolumn{11}{l}{\textit{Panel A: Depression}} \\"

eststo clear
local col 1	
	
	*Loop over countries
	levels d_country, local(d_country)
	local forlab: value label d_country
	foreach country of local d_country{
	local country_label: label `forlab' `country'
		
	eststo col`col++': reg depression i.d_unemp d_unemp#c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	if (c_lms_new==1 | c_lms_new==2 | c_lms_new==3) /// 
	& ragender==`x' & ragey>=50 & ragey<=70 & d_country==`country', vce(cluster nmergeid)
	
	local obs = e(N)
	estadd scalar obs= `obs'
	
	}
	
	esttab using $out_files/`file', append b(a2) se(a2) keep(0.d_unemp#c.ragey 1.d_unemp#c.ragey 1.d_unemp _cons) label nolines nogaps nomtitles compress fragment nonumbers mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) stats(obs, labels("Observations")) coeflabels(0.d_unemp#c.ragey "Employed/Retired $\times$ Age" 1.d_unemp#c.ragey "Unemployed $\times$ Age" 1.d_unemp "Unemployed")
	
	writeln $out_files/`file' "\midrule"

*EURO-D
writeln $out_files/`file' "\multicolumn{11}{l}{\textit{Panel B: EURO-D score}} \\"

eststo clear
local col 1	
	
	*Loop over countries
	levels d_country, local(d_country)
	local forlab: value label d_country
	foreach country of local d_country{
	local country_label: label `forlab' `country'
					
	eststo col`col++': reg reurod i.d_unemp d_unemp#c.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	if (c_lms_new==1 | c_lms_new==2 | c_lms_new==3) /// 
	& ragender==`x' & ragey>=50 & ragey<=70 & d_country==`country', vce(cluster nmergeid)
	
	local obs = e(N)
	estadd scalar obs= `obs'
	
	}
	
	esttab using $out_files/`file', append b(a2) se(a2) keep(0.d_unemp#c.ragey 1.d_unemp#c.ragey 1.d_unemp _cons) label nolines nogaps nomtitles compress fragment nonumbers mlabels(none) starlevels(* 0.10 ** 0.05 *** 0.01) stats(obs, labels("Observations")) coeflabels(0.d_unemp#c.ragey "Employed/Retired $\times$ Age" 1.d_unemp#c.ragey "Unemployed $\times$ Age" 1.d_unemp "Unemployed")
	
	writeln $out_files/`file' "\midrule"
	
}

********************************************************************************
* Appendix Table F.2: Summary statistics by labor-market status for men aged 50-70 
********************************************************************************
clear all
use $working_data/SHARE_JOLE_replication.dta, clear

local file summary_table.tex

local binary "depression","educ_no","educ_primary","educ_sec_low","educ_sec_up", ///
	"educ_post_sec","educ_tert_first","educ_tert_sec", "d_early_ret"
	local binarytwo "d_normal_ret","mar_married","mar_partner","mar_separ","mar_divor","mar_widow","mar_never"
	local outcomes  depression reurod
	local controls ragey  hittot radla riadlza
	local education educ_no educ_primary educ_sec_low ///
	educ_sec_up educ_post_sec educ_tert_first educ_tert_sec
	local maritalstat mar_married mar_partner mar_separ mar_divor mar_widow mar_never
	
	g full  = 1 & ragender==1 
	g smpl1 = full & c_lms_new==1 & ragey>=50 & ragey<=70
	g smpl2 = full & c_lms_new==2 & ragey>=50 & ragey<=70
	g smpl3 = full & c_lms_new==3 & ragey>=50 & ragey<=70
	g smpl4 = full & c_lms_new==4 & ragey>=50 & ragey<=70
	g smpl5 = full & c_lms_new==9 & ragey>=50 & ragey<=70
	
	local conditions "smpl1" "smpl2" "smpl3" "smpl4" "smpl5" 
	
	local colhead1 & Men 
	local colhead2 & Retired & Employed & Unemployed & Disabled & Other

	local cols : word count "`conditions'"
	local cols = `cols'+1

	writeln $out_files/`file' "\begin{table}[H]\centering "
	writeln $out_files/`file' "\caption{Pooled sample - Summary statistics for the pooled sample \label{tab:summary_table}}"
	writeln $out_files/`file' "\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}"
	writeln $out_files/`file' "\begin{threeparttable} "
	writeln $out_files/`file' "\begin{tabular}{l *{`=`cols'-1'}{c}} "		
	writeln $out_files/`file' "\toprule  "
	writeln $out_files/`file' " `colhead1'    \\"	
	writeln $out_files/`file' " `colhead2'    \\"	
	local paneltitles "Outcome "Characteristics" "Education (ISCED-97)" "Marital Status""
	local panels outcomes controls education maritalstat
	foreach panel in `panels'  {
		local ppos: list posof "`panel'" in local(panels)
		local paneltit: word `ppos' of `paneltitles'
		writeln $out_files/`file' "\midrule"		
		writeln $out_files/`file' "\multicolumn{`cols'}{l}{\textbf{`size'{`paneltit'}}} \\ "
		foreach v in ``panel'' {
			local vl : variable label `v'
			local linemean `vl' 
			local linesd
			foreach c in "`conditions'" {
				if "`c'"=="full" & (("`panel'"=="dwscontrols") | ("`=substr("`v'",1,2)'"=="lj")) {
					local linemean  `linemean' & 
					local linesd    `linesd'   & 
					continue
				}
				qui sum `v' if `c'
				sigdigits a2 `=r(mean)'
				local mean: disp `fmt' `=r(mean)'
				local linemean  `linemean' & `mean'
				sigdigits a2 `=r(sd)'
				local sd: disp `fmt' `=r(sd)'
				local linesd    `linesd'   &  `=cond(`sd'!=.,"[`=ltrim("`sd'")']","`sd'")' //  ,`=r(N)'
			}
			writeln $out_files/`file' "`linemean' \\"
			if !inlist("`v'","`binary'") & !inlist("`v'","`binarytwo'") writeln $out_files/`file' "`linesd'   \\"
			
		}	
	}
	writeln $out_files/`file' "\midrule  "	
	local j 1
	local countline 
	foreach c in "`conditions'" {
		count if `c'
		local count`j' = r(N)
		local countline `countline' & `count`j++''
	}
	writeln $out_files/`file' "Number of observations `countline' \\ "	

	writeln $out_files/`file' "\bottomrule "
	writeln $out_files/`file' "\end{tabular}"
	writeln $out_files/`file' "\end{threeparttable} "
	writeln $out_files/`file' "\end{table}"

********************************************************************************
* Appendix Tables F.3 and F.4: Cross-tabulation of the number of observations (N) by self-reported and income-based labor-market status between ages 50 - 70 
********************************************************************************
clear all
use $working_data/SHARE_JOLE_replication.dta, clear

*Men (Appendix Table F.3)
latab self_lf c_lms_new if ragender==1 & ragey>=50 & ragey<=70, tf($out_files/lms_cross_men) replace

*Women (Appendix Table F.4)
latab self_lf c_lms_new if ragender==2 & ragey>=50 & ragey<=70, tf($out_files/lms_cross_women) replace
	

********************************************************************************
* Appendix Table F.5: Correlates of mental well-being 
********************************************************************************
clear all
use $working_data/SHARE_JOLE_replication.dta, clear

local file ushape_correlates.tex

eststo clear
local col 1	

	*Depression 
	eststo col`col++': reg depression i.ragender i.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	ib2.c_lms_new ///
	if ragey>=50 & ragey<=70 , vce(cluster nmergeid)
	
	
	*EURO-D
	eststo col`col++': reg reurod i.ragender i.ragey ///
	hh_inc i.raedisced i.rmstat i.wave ///
	i.isocountry ///
	ib2.c_lms_new ///
	if ragey>=50 & ragey<=70 , vce(cluster nmergeid)
	
	esttab using $out_files/`file', drop(*.isocountry *.wave *.ragey) ///
					nogaps booktabs replace se label ar2 fragment nonumbers ///
					b(a3) star(* 0.10 ** 0.05 *** 0.01) ///
					title("Self-reported retirement") ///
					coeflabels(2.ragender "Gender (female=1)" ///
					1.d_unemp#1.d_rep_retired "Unemployed $\times$ Reported ret." ///
					d_early "Eligibility" 1.d_unemp#1.d_early "Unemployed $\times$ Eligibility" ///
					1.raedisced "Primary" 2.raedisced "Lower secondary" 3.raedisced "Upper secondary" ///
					4.raedisced "Post-secondary" 5.raedisced "First stage tertiary" 6.raedisced "Second stage tertiary" ///
					1.rmstat "Married" 3.rmstat "Partnered" 4.rmstat "Separated" 5.rmstat "Divorced" ///
					7.rmstat "Widowed" 8.rmstat "Never married" ///
					hh_inc "HH income ('0 000 EUR)" log_hhinc "log(HH income)" radla "ADL" riadlza "IADL" ///
					_cons "Constant") /// 
					nomtitles substitute(\_ _) 
				
********************************************************************************
* Appendix Figure G.1: Fraction of retirees by country and survey wave over ages 50-70 for men 
********************************************************************************	
*Note: below code also generates figures for women, which are not included in the paper
				
clear all
use $working_data/SHARE_JOLE_replication.dta, clear
					
*Create retirement fraction variable
levels ragender, local(ragender)
local forgender: value label ragender
*Loop over gender
foreach x of local ragender{
local gender_label: label `forgender' `x'
	
	levels wave, local(wave)
	foreach wave in `wave'{	
		
		*Create matrix with fraction of retirees
		matrix Ret_`gender_label'_`wave' = J(31,11,-99)
		
			*Loop over countries
			levels d_country, local(d_country)
			foreach country in `d_country'{
				forvalues i=50(1)80{
				
				local age_condition="inrange(ragey, `i'-1, `i'+1)" 

				matrix Ret_`gender_label'_`wave'[`i'-49,1]=`i'
				sum c_lms_new if ragender==`x' & `age_condition' & d_country==`country' & wave==`wave' & (c_lms_new==1 | c_lms_new==2)
				scalar total=r(N)
				sum c_lms_new if ragender==`x' & `age_condition' & d_country==`country' & c_lms_new==1 & wave==`wave'
				scalar retired=r(N)
				matrix Ret_`gender_label'_`wave'[`i'-49,`country'+1]=retired/total
				}
			}

			matrix list Ret_`gender_label'_`wave'
			
			matrix colnames Ret_`gender_label'_`wave' = "Age" "Austria" "Belgium" "Denmark" "France" ///
			"Germany" "Italy" "Netherlands" "Spain" "Sweden" "Switzerland"
			
			svmat Ret_`gender_label'_`wave', names(matcol)	
		
	}
	
	*Create graphs
	levels d_country, local(d_country)
	local forlab: value label d_country
		foreach country of local d_country{
		local country_label: label `forlab' `country'
		
				*Select early retirement ages for plots (men only)
				if "`country_label'"=="Austria" & "`gender_label'"=="Men" local early_line="xline(60, lcolor(gray*0.40) lpattern(solid))" 
				if "`country_label'"=="Belgium" & "`gender_label'"=="Men" local early_line="xline(60(.01)63, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="Denmark" & "`gender_label'"=="Men" local early_line="xline(60(.01)64, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="France" & "`gender_label'"=="Men" local early_line="xline(56, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="Germany" & "`gender_label'"=="Men" local early_line="xline(60(.01)63, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="Italy" & "`gender_label'"=="Men" local early_line="xline(57(.01)60, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="Netherlands" & "`gender_label'"=="Men" local early_line="xline(60(.01)62, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="Spain" & "`gender_label'"=="Men" local early_line="xline(60(.01)63, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="Sweden" & "`gender_label'"=="Men" local early_line="xline(61, lcolor(gray*0.40) lpattern(solid))"
				if "`country_label'"=="Switzerland" & "`gender_label'"=="Men" local early_line="xline(63, lcolor(gray*0.40) lpattern(solid))"	
		
		line Ret_`gender_label'_1`country_label' Ret_`gender_label'_1Age if Ret_`gender_label'_1Age<=70, lcolor(black) lpattern(solid) || ///
		line Ret_`gender_label'_2`country_label' Ret_`gender_label'_2Age if Ret_`gender_label'_1Age<=70, lcolor(black) lpattern(dash) || ///
		line Ret_`gender_label'_4`country_label' Ret_`gender_label'_4Age if Ret_`gender_label'_1Age<=70, lcolor(black) lpattern(shortdash) || ///
		line Ret_`gender_label'_5`country_label' Ret_`gender_label'_5Age if Ret_`gender_label'_1Age<=70, lcolor(gs10) lpattern(solid) || ///
		line Ret_`gender_label'_6`country_label' Ret_`gender_label'_6Age if Ret_`gender_label'_1Age<=70, lcolor(gs10) lpattern(dash) ///
		xtitle("Age") xlabel(50(2)70) `early_line' ///
		ytitle("Fraction of retirees", axis(1)) ylabel(0(.2)1, angle(0) axis(1)) ///
		legend(label(1 "Wave 1 (2004/5)") label(2 "Wave 2 (2006/7)") label(3 "Wave 4 (2011/12)") label(4 "Wave 5 (2013)") label(5 "Wave 6 (2015)") rows(2)) ///
		legend(off) graphregion(color(white)) ///
		title(`country_label', color(black)) ///
		name(Ret_`country_label'_`gender_label')
	
	}
	
	*Combine graphs 
	grc1leg Ret_Austria_`gender_label' Ret_Belgium_`gender_label' ///
	Ret_Denmark_`gender_label' Ret_France_`gender_label' ///
	Ret_Germany_`gender_label' ///
	Ret_Italy_`gender_label' Ret_Netherlands_`gender_label' Ret_Spain_`gender_label' ///
	Ret_Sweden_`gender_label' Ret_Switzerland_`gender_label' ///
	, legendfrom(Ret_Austria_`gender_label') ///
	row(4) xsize(9) graphregion(color(white)) name(RetFrac_`gender_label'_paper) 
				
	graph display RetFrac_`gender_label'_paper,  xsize(3.5) ysize(3.8) scale(0.85)
	
	graph export "$out_files/RetFrac_`gender_label'_paper.pdf", replace	
	}
	
	
